1-D Heat Equation#
강좌: 기초 전산유체역학
Heat Equation (1-D)#
열전도에 의한 열유속 (Heat Flux)는 Fourier Law에 의해 표현할 수 있다.
여기서 \(T\)는 온도, \(k\) 는 열전도 계수이고, \(\dot{q}\) 는 열유속 (\(W/m^2\)) 이다.
단위 길이 \(\Delta x\) 에 대해서 내부에너지 변화는 열유속의 구배와 같다.
이를 정리하면 1차원 Heat equation 이다.
Fig. 10 Heat Conduction (image from wikimedia.org)#
Finite Difference Method#
Linear Wave 방정식 해석과 같이 Finite Difference 를 이용해서 이 편미분 방정식을 차분한다.
FTCS (Forward difference in Time, Central difference in Space)#
여기서 시간과 공간을 각각 M, N 개로 나누고, \((n, i)\) 점에서 값을 다음과 같이 간단히 표시하자

Fig. 11 FTCS Stencil#
그 결과 위 차분식 (FDE, Finite Difference Equation)은 다음과 같이 표현할 수 있다.
계산 격자 구성, Solution array 구성#
계산 영역을 \(n_x + 1\)개의 점으로 나누어보자.
즉 격자점은 다음과 같다.

Fig. 12 Grid#
첫번째 격자점과 마지막 격자점은 경계조건으로 부여할 수 있다.
경계 조건을 위해서 Solution array는 경계조건을 포함해서 격자점 개수와 같이 \(n_x+1\) 으로 구성한다.
계산은 양 끝점을 제외한 \(T_1, T_2,... T_{n_x}\) 에 대해서 수행한다.
예제#
\(\alpha=0.1\) 이고 \(x\in [0, 1]\) 에 대해서 해석한다.
초기 조건 : \(T(x,0) = 100\)
경계 조건 : \(T(0, t) = 100\), \(T(1, t) = 300\)
\(\Delta t =0.01, \Delta x=0.1\) 에 \(t=1\) 일때 온도 분포를 구하시오.
참고 Exact Solution은 변수 분리법에 따라 다음과 같다.
%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
plt.style.use('ggplot')
plt.rcParams['figure.dpi'] = 150
def cs(nx, u, du, adtdx2):
"""
Central difference in space
Parameters
----------
nx : integer
number of elements in solution array
u : array
solution
du : array
difference of current and next solution
adtdx2 : float
constant
"""
# DIY
# du_i = a*dt/dx^2*(u_{i-1} - 2*u_i + u_{i+1}) for i in 1, 2..., nx-2
# Conditions
alpha = 0.1
ti = 100
ts = 300
t_target = 1.0
dt = 0.01
# Make grid
nx = 11
x = np.linspace(0, 1, nx)
dx = np.diff(x)[0]
# Const
adx2 = alpha / dx**2
# Solution array
u = np.ones_like(x)*ti
du = np.zeros_like(u)
# Calculation
t = 0
while abs(t - t_target) > 1e-10:
# Adjust time step to reach target time
dt = min(dt, t_target - t)
# Apply bc
u[0] = ti
u[-1] = ts
# Spatial discretization
cs(nx, u, du, adx2*dt)
# Update solution and time
u += du
t += dt
# Exact solution
dts = ts - ti
u_exact = ti + dts*x + np.sum([
2*dts*(-1)**n / (n*np.pi) *np.exp(-n**2*np.pi**2*alpha*t)*np.sin(n*np.pi*x)
for n in range(1, 10)
], axis=0)
# Visualization
plt.plot(x, u, marker='x')
plt.plot(x, u_exact)
plt.legend(['Computation', 'Exact'])
# Compute Error
err = np.linalg.norm(u[1:-1] - u_exact[1:-1])
err = u - u_exact
err = np.sqrt(np.sum(err**2) / nx)
print('Error: {:.5f}'.format(err))
Error: 0.16054

Accuracy, Stability Convergence#
정확도#
FTCS 의 Truncation Error 는 \(O((\Delta t), (\Delta x)^2)\) 이다. 즉 시간 차분에 대해서는 1차 정확도, 공간 차분에 대해서 2차 정확도를 갖는다.
Exact solution 과의 오차는 다음과 같이 계산한다.
선형 편미분 방정식이므로, Linear Wave 방정식과 같이 Consistency, stability 분석이 가능하다.
또한 이를 통해 수치 해석 결과가 이론 해로 수렴 (Convergence) 하는지 확인할 수 있다.
von Neumann Stability#
FDE (Finite Difference Equation) 의 완전해를 \(D\) 라 했을 때 수치해 \(N\)은 다음과 같이 나타낼 수 있다.
여기서 \(\epsilon\) 은 round-off 에 의한 error 이다.
이를 차분식에 적용하면, \(D\) 은 차분식을 만족하므로 사라지고, \(\epsilon\) 만 남는다.
이 오차를 다음과 같은 형태로 표현하자.
이를 차분식에 적용하면
정리하면
\(|\sigma| < 1\) 을 만족하기 위해서는
즉
Convergence#
Consistency와 Stability를 만족하므로 Lax Equivalance Theorem에 의해 수치해가 이론해에 수렴함을 알 수 있다.
Du Fort-Frankel Method#
FTCS에서 시간에 대해서도 Central Difference를 하면
이 기법은 Unconditionally unstable 하다.
수치 안정성을 확보하기 위해 우변 항을 변화하면 Du Fort-Frankel 기법이다.
이 기법은 Unconditionally stable 하다.
Talyor expansion을 이용해 Truncation Error을 확인하면
\(\Delta t \rightarrow 0\), \(\Delta x \rightarrow 0\) 으로 점근하더라도 Truncation Error는 사라지지 않는다.
즉 Consistency를 만족하지 않는다.
실습#
FTCS 기법을 완성하시오.
FTCS 기법에 대해서 \(\Delta t\) 를 바꿔가면서 von Neumann Stability Analysis 결과를 확인하시오.
\(\Delta t=10^{-5}\) 로 매우 작게 정하자. 이때 \(\Delta x=0.1, 0.05, 0.025, 0.0125\) 로 바꿔가면서 오차를 구하시오. 공간에 대해 2차 정확도를 만족하는지 확인하시오.
\(\Delta x=0.05\) (\(n_x=21\)) 로 정하자. 이때 \(\Delta t=0.01, 0.005, 0.0025, 0.00125\) 로 바꿔가면서 오차를 구하시오. 시간에 대해 1차 정확도를 만족하는지 확인하시오.